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Abstract 

Background: During euloryotic genome evolution, tandem gene duplication is the most frequent event giving rise to 
clustered gene families. However, how expression divergence between tandemly duplicated genes has emerged and 
maintained remain unclear. In particular, it is unknown if epigenetic regulators have been involved in the process. 

Results: We demonstrate that CCCTC-binding factor (CTCF), the master epigenetic regulator and the only l<nown 
insulator protein in humans, has played a predominant role in generating divergence in both expression profiles 
and expression levels between adjacent paralogs in the human genome. This phenomenon was not observed for 
non-paralogous adjacent genes. After tandem duplication events, CTCF-binding sites gradually accumulate 
between paralogs. This trend was more prominent for genes involved in particular functions. 

Conclusions: The accumulation of CTCF-binding sites drives expression divergence of tandemly duplicated genes. 
This process is likely targeted by natural selection. Our study reveals the importance of CTCF to the evolution of 
animal diversity and complexity. 




Genomics 



Background 

Gene duplication is a major driver for the emergence of 
organismal complexity and evolutionary innovations 
[1,2]. Tandem duplication is the most common route to 
the formation of clustered paralogous genes [3,4]. A 
newly duplicated gene must diverge from its progenitor 
gene in coding sequence or expression, or it will degen- 
erate due to redundancy [1,5,6]. To acquire novel tran- 
scription patterns, tandemly duplicated genes need to 
interrupt expression similarity due to shared upstream 
cw-elements upon origin [3,7] and transcriptional inter- 
ference contributed by physical proximity [8]. Despite 
the challenges, functionally important gene clusters con- 
sisting of paralogs with distinct expression patterns are 
found in a wide range of species, including human 
[9,10]. Therefore, it is important to understand the 



origin and maintenance of expression divergence 
between tandem paralogs. 

CCCTC-binding factor (CTCF), the only known 
human insulator protein, plays a master role in deter- 
mining the transcriptional landscape of genomes. When 
bound to insulator sequences (CTCF-binding sites), 
CTCF can prevent repressive heterochromatin from 
spreading into neighboring regions [11,12]. In addition, 
CTCF interferes with enhancer-promoter communica- 
tion [13] and guides long-range chromatin interactions 
[14]. Because changes in epigenetic marks play roles in 
regulatory divergence of duplicated genes [15-17], we 
hypothesize that CTCF also plays a role. Using human 
RNA-seq and ChlP-seq data (see Methods), we exam- 
ined how CTCF drives regulatory divergence of dupli- 
cated genes, especially tandemly arrayed paralogs. 
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Results and discussion 

In the genome, adjacent genes are coexpressed due to a 
common origin [18], cofunctionality [19], or deleterious 
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transcriptional interference [8]. Using RNA-seq data 
from six human adult tissues [20], we measured gene 
expression dissimilarity between adjacent genes in terms 
of expression profile or expression level using ExpDi_^ or 
ExpD^ac (see Methods), respectively. ExpDi.^ focused on 
changes in the shape across the tissue dimension, while 
ExpD^^c focused on summed changes in abundance. We 
studied divergences in both profiles and abundances 
because previous studies have shown that the observa- 
tions from the two aspects did not necessarily produce 
consistent results [21], possibly due to different underly- 
ing mechanisms controlling the properties [22]. 

Mammalian genes close to each other have similar 
expression profiles [8,19]. CTCF-binding sites can pre- 
vent undesirable crosstalk between active and inactive 
genomic regions [23]. DNA methylation upstream of a 
gene can inhibit its transcriptional initiation [24]. To 
determine if the effect of intergenic distance, the num- 
ber of CTCF-binding sites, and DNA methylation is 
related to shared evolutionary origin, we measured 
ExpDi^r and ExpD^ac in non-paralogous adjacent genes. 
Non-paralogous gene pairs with a longer intergenic dis- 
tance {d, calculated as the distance in nucleotides 
between the transcriptional start sites, TSS), more over- 
lapping CTCF-binding sites {iCTCF, see Methods), or a 
larger difference in upstream DNA methylation 
(AQjGo/e. see Methods) had greater ExpDi^^ or ExpD^^^ 
(see p in Table 1). Although d, #CTCF, and ACpGo/E 
were interrelated (d vs. #CTCF: p =0.521, P<10'^°''; 
#CTCFvs. ACpGo/E- p=0.075, P<10'*°; AQjGq/e vs. d: 
p=0.166, P<10'^^^), partial correlation analyses (see 
Methods) suggested that the effects of d, #CTCF, and 
ACpGo/E on generating expression dissimilarity between 
adjacent genes was independent, d had the largest direct 
influence on ExpD-i-v:, and ACpGo/E had the largest 
direct influence on ExpD^^^. #CTCF had a weak but sig- 
nificant effect on ExpDx-r and an intermediate effect on 
ExpD^^c (see pp in Table 1). Although CTCF binding 
can vary among cell types [25], when we define #CTCF 
using joint CTCF ChlP-seq peaks (see Methods) instead 
of overlapping CTCF peaks, the results did not change 
(Table SI in Additional file 1). 



Mechanisms generating coding sequence divergence, 
such as change in protein structure or splicing, between 
paralogous genes have been intensively investigated 
[2,26-28]. By contrast, mechanisms generating expres- 
sion divergence have garnered less attention. Although 
#CTCF was not the strongest determinant of expression 
divergence in non-paralogous adjacent gene pairs 
(Table 1), for adjacent paralogs, #CTCF had the greatest 
direct influence on ExpDi_^ and ExpD^^^i followed by 
ACpGo/E (see Pp in Table 2). Repeating analysis based on 
an independently generated and unpublished RNA-seq 
dataset including 16 human tissues (lUumina BodyMap 
2.0 project, see Methods) produced a result consistent 
with Table 1 and 2 (Table S2 in Additional file 1), 
suggesting the robustness of the pattern found. 

We classified adjacent paralogs into three groups 
according the orientations: head-to-head, head-to-tail 
and tail-to-tail [29]. Although reduced sample sizes 
resulted in reduced statistical significance, the pattern of 
the greatest impact #CTCF on expression divergences 
held regardless of the orientation of paralogs (Table 3). 
Using microarray expression data, a previous study 
found that intergenic distance was related to expression 
divergence for tandemly arrayed paralogs [30]. However, 
microarray data is known to have cross-hybridization 
related biases [17], and we found no direct association 
between intergenic distance and expression divergence 
for tandemly arrayed paralogs after controlling for 
#CTCF and ACpGo/E (see pp of d vs. ExpDi^,- or 
ExpDEuc, Table 2). Taken together, these results suggest 
that CTCF-binding sites play a very significant, if not 
primary, role in driving expression divergence of tan- 
demly duplicated genes. 

There are two hypotheses to explain the influence of 
CTCF-binding sites in driving expression divergence of 
tandem paralogs. First, tandem paralogs that arose in 
genomic regions with high densities of CTCF-binding 
sites nearby are more likely to be preserved due to 
immediate independence of gene regulation. Second, 
CTCF-binding sites accumulated over time to enhance 
independent gene regulation of tandem paralogs, espe- 
cially those have been functionally diverged. If the first 



Table 1 Rank correlations (p) and partial rank correlations (pp) of the genomic properties with expression 
dissimilarities (measured by ExpD^.r or ExpDEuc) of the non-paralogous adjacent genes 



Genomic properties ' 




B<PDl-r 




BfpDEuc 


p (P-value) 


Pp (P-value) 


p (P-value) *■ 


Pp (P-value) 


d 


0.180 (<10"^^') 


0.130 (<10""^) 


0.112 (<10"^^) 


0.040 (<10"") 


tCTCF 


0.120 (<10"''^) 


0.032 (<10"') 


0.104 (<10"") 


0.058 (<10"^^) 




0.0/ ,' i vO 




O.:o6 ;-. iO-"'"i 


0,182 0 0 ■'-'■') 



^ "d", intergenic distance; "#CJCF", number of overlapping CTCF-binding sites; "ACpGo/E - difference in upstream DNA methylation. 
^ P values show the probabilities of the observations under the hypothesis of no correlation. 

^ Spearman's partial correlation coefficient pp is computed by controlling for the other two genomic properties listed in ^ 



Liao and Chang BMC Genomics 2014, 15(Suppl 1):S8 
httpy/www.biomedcentral.coni/1 471 -21 64/1 5/S1 /S8 



Page 3 of 7 



Table 2 Rank correlations (p) and partial rank correlations (pp) of the examined genomic properties with expression 
dissimilarities (measured by £xpD,.r or ExpD^uc) of the paralogous adjacent genes 



Genomic properties ' 




ExpD^., 




ExpDeuc 


p (P-vaiue) 


Pp (P-value) 


p (P-value) 


Pp (P-vaiue) 


d 


0.039 (0.164) 


-0.011 (0.707) 


0.100 {<10"^) 


0.021 (0.461) 


tacF 


0.111 (<10"^) 


0.103 (<10"') 


0.187 (<10"^°) 


0.159 (<10"^) 


ACpGo/E 


0.069 (0.015) 


0.054 (0.023) 


0.102 (<10"^) 


0.094 (<10"^) 



" "d", intergenic distance; "tCTCF", number of overlapping CTCF-binding sites; "MlpGo/i', difference in upstream DNA methylation, 
P values show the probabilities of the observations under the hypothesis of no correlation. 

Spearman's partial correlation coefficient pp is computed by controlling for the other two genomic properties listed in ' 



hypothesis is correct, there would be no correlation 
between the divergence time of paralogs (measured by 
ds or Tphy, see Methods) and #CTCF. However, we 
observed a positive correlation between ds or Tp^y and 
#CTCF (Table 4), suggesting that the second hypothesis 
is correct. #CTCF had stronger rank correlation (p), 
which is positive, to d^ (or Tphy) than d and AQ^Gq/e 
(Table 4). Partial correlation analyses further suggested 
that the increase of #CTCF with d^, (or Tphy) is not 
caused by the change in d or AQjGq/e or their com- 
bined effect over time (see pp in Table 4). By contrast, 
the insignificant partial correlation of d with respect to 
ds (or Tphy) indicated that the increase d over time can 
be explained by the increase in the number of CTCF- 
binding sites and the associated changes in DNA 
methylation. 

To determine whether tandem paralogs with specific 
functions tend to have a greater number of intervening 
CTCF-binding sites, we performed enrichment analyses 
in Gene Ontology (GO) terms. To eliminate the poten- 
tial effect of duplicability in GO analysis [4], paralogous 
and non-paralogous gene pairs were analyzed separately. 



To control for the potential effect of gene density [31], 
#CTCF/d was used instead of #CTCF (although using 
#CTCF produced a consistent result, which is not 
shown). Gene pairs in the top quartile of CTCF-binding 
site density {#CTCF/d) were compared against the bot- 
tom three quartiles. Enriched/depleted GO terms for the 
paralog group were substantially different from those in 
the non-paralog group (Table 5). Only two GO terms 
(GO:0010033, response to organic substance; 
GO:0031012, extracellular matrix) had the same enrich- 
ment status ("enriched") in both groups. Tandemly 
duplicated genes with a higher density of intervening 
CTCF-binding sites tended to specifically encode pro- 
teins involved in gene expression (GO:0010467, 
GO:0008134), metabolic processes (GO:0019222, 
GO:0006139, GO:0050790), or cellular processes 
(GO:0050794, GO:0051128, GO:0044249) through DNA 
binding (GO:0003677), SMAD binding (GO:0046332), 
growth factor binding (GO:0019955), or kinase interac- 
tion (GO:0019210, GO:0019887) in receptor complexes 
(GO:0043235) or intracellular regions (GO:0031012) 
(Table 5). This result implied that the densities of CTCF 



Table 3 Rank correlations (p) and partial rank correlations (pp) of the examined genomic properties with expression 
dissimilarities (measured by ExpO^., or ExpD^„^ of the adjacent paralogous genes of different orientations. 



Genomic properties ' ExpO^., ExpD^uc 





p (P-value) 


Pp (P-value) "'^ 


p (P-value) 


Pp (P-value) 


Head-to-Head 










d 


0.1 65 (0.072) 


0.106 (0.250) 


0.150 (0.109) 


0.036 (0.700) 


ttacF 


0.213 (0.020) 


0.173 (0.058) 


0352 (<10''*) 


0.327 (<10'^) 


ACpGo/E 


-0.024 (0.794) 


0.001 (0.991) 


-0.071 (0.443) 


-0.054 (0.561) 


Head-to-Tail 










d 


0.005 (0.870) 


-0.050 (0.110) 


0.089 (<10"^) 


0.016 (0.607) 


iCTCF 


0.100 (<10"^) 


0.110 (<10"^) 


0.151 (<10"=) 


0.127 (<10"') 


ACpGo/E 


0.083 (<10'^) 


0.085 (<10"^) 


0.109 (<10"') 


0.104 (<10"^) 


Tail-to-Tail 










d 


0.139 (0.110) 


0.105 (0.230) 


0.145 (0.096) 


0.046 (0.602) 


ttCTCF 


0.113 (0.195) 


0.061 (0.482) 


0.271 (<10"^) 


0.202 (0.019) 


ACp6o/E 


0.050 (0.569) 


0.018 (0.0835) 


0.212 (0.014) 


0.153 (0.078) 



^ "d", intergenic distance; "#CrCF", number of overlapping CTCF-binding sites; "ACpGo/E, difference in upstream DNA methylation. 
^ P values show the probabilities of the observations under the hypothesis of no correlation. 

Spearman's partial correlation coefficient pp is computed by controlling for the other two genomic properties listed in ^ 
^ There were 120, 1003 and 133 head-to-head, head-to-tail and tail-to-tail adjacent paralogs, respectively. 
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Table 4 Rank correlations (p) and partial rank correlations (pp) of the examined genomic properties with the 
divergence time of the paralogous adjacent genes. 

Genomic properties ' ds Tphy 

p (P-value) Pp (P-value) p (P-value) pp (P-value) ''■^ 

All adjacent paralogs 

d 0.138 (<10"^) 0.057 (0.056) 0.141 (<10"**) 0.053(0.060) 

^CTCF 0.225 (<10"^^) 0.189 (<10''') 0.219 (<10"'*) 0.178 (< 10"^) 

ACp6o/E 0.097 {< 10"^) 0.089 (< 10'^) 0.084 (< 10'^) 0.073 (< 10"^) 
Adjacent paralogs associated with GO terms in which high #CTCF/d genes were specifically enriched 

d 0.086 (0.152) -0.070 (0.245) 0.108(0.072) -0.038 (0.526) 

#aCF 0.292 {<10"'^) 0.276 (<10'^) 0.273 (< 10'^) 0.243 (< 10"*) 

ACpGo/E 0.122(0.043) 0.079 (0.186) 0,132 (0.027) 0,092(0,125) 

" "d", intergenic distance; "ItCTCF", number of overlapping CTCF-binding sites; "ACpGo/E, difference in upstream DNA methylation. 
^ P values show the probabilities of the observations under the hypothesis of no correlation. 
Spearman's partial correlation coefficient pp is computed by controlling for the other two genomic properties listed in ' 



binding sites between tandem paralogs were not con- 
tributed by the genomic background. 

Genomic regions of high densities of CTCF binding 
sites can emerge through stochastic evolutionary pro- 
cesses. To examine if high CTCF-binding site density 
between tandem duplicates is the outcome of gradual 
CTCF-binding site accumulation by natural selection, 
we focused on the subset of 278 paralogous pairs in 
which at least one paralog had one of the abovemen- 
tioned enriched GO categories. In this subset of tandem 
paralogs, Pp of #CTCF vs. ds (or Tp^y) after controlling 
for d and ACpGo/E was stronger {#CTCF vs. ds: 
Pp=0.276, P<10'^; #CTCFvs. T^^y. pp=0.243, P<10'*) 
than that observed in the full set of adjacent paralogs 
{#CTCF vs. ds. pp=0.189, P<10'^; #CrCF vs. T^^y. 
Pp=0.178, P<10''') (Table 4). When #CTCF was defined 
by joint CTCF ChlP-seq peaks rather than overlapping 
CTCF peaks, the results were similar (Table S3 in Addi- 
tional file 1). Therefore, the trend to accumulate CTCF- 
binding sites following a tandem duplication event 
resulted in more CTCF-binding sites (both absolute 
number and density) between paralogs, especially for 
those with enriched GO-categories shown in Table 5. 

Conclusions 

Combining human genomic and transcriptomic data, 
this study demonstrates that CTCF and its binding sites 
play a major role driving the expression evolution of 
tandemly duplicated genes. Following tandem duplica- 
tion events, CTCF-binding sites gradually accumulate 
between the paralogs to increase their divergences in 
expression profile and their divergences in expression 
level. The role of CTCF-binding sites is not limited to 
the insulation of DNA methylation domains, because 
the effects of #CTCF on ExpDi.^ (or ExpO^^^) were still 
significant even after controlling for ACpGo/E (Table 2). 
Thus, CTCF, a conserved regulatory protein [32], affects 



the expression evolution of adjacent genes in genomes 
from flies [33] to humans and is important for the evo- 
lution of organismal complexity in animals. 

Methods 

Annotation of the human genome (Ensembl v72), includ- 
ing gene coordinates, TSS, and paralog divergence, was 
retrieved through BioMart (http://www.biomart.org/) 
[34]. Orientation for each paralogous pair (head- to-head, 
head-to-tail, or tail-to-tail) was determined using strand 
information. For each paralogous pair, the rate of synon- 
ymous changes (ds) was calculated using PAML [35], and 
the phylogenetic age (Tphy) was assigned according to 
Table S4 in Additional file 1 based on Ensembl's annota- 
tion of the most recent common ancestor [36]. A smaller 
ds or Tphy indicated a more recent divergence time. The 
RNA-seq-based gene expression signals in the human 
brain, cerebellum, heart, kidney, liver and testis 
(GSE30352) [20] were obtained from NCBI Gene Expres- 
sion Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo/ 
). The raw reads of RNA-seq data in 16 human tissues 
(adipose, adrenal, brain, breast, colon, heart, kidney, liver, 
lung, lymph node, ovary, prostate, skeletal muscle, testes, 
thyroid, white blood cells) by Illumina BodyMap 2.0 pro- 
ject were downloaded from GSE30611 of GEO and were 
processed following our previous studies [17,37] to 
obtain expression signals. Upstream DNA methylation 
for a gene was calculated based on the 500 nucleotides 
upstream of the TSS by CpGo/E=J'cpG /{Pc x Pg)> where 
PcpG' Pc> and Pq are the frequencies of CpG dinucleo- 
tides, C nucleotides, and G nucleotides, respectively 
[17,38]. The difference in upstream DNA methylation, 
ACpGo/E, was calculated by the absolute value of the dif- 
ference in CpGo/E between the two adjacent genes com- 
pared. In total, we obtained 30,164 non-paralogous and 
1,256 paralogous gene pairs based on 32,164 human 
genes with detectable expression and estimable CpGo/E- 
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Table 5 Enriched/depleted GO terms at level 4 for genes with high #CTCF/d flanking regions 



Functional 
categories 



Duplicate with 25% top #CTCF/d vs. rest of duplicates 



Non-duplicates with 25% top #CTCF/d vs. 
rest of non-duplicates 



GO terms ' (t or P-value ^) 



GO terms (t or P-value 



iVIolecular functions 

DNA binding 

(GO;0003677) (t, P<0.05) 
transcription factor binding 

(GO:0008134) (f, P<10"^) 

SiVIAD binding 

(GO:0046332) (j, P<0.05) 

growth factor binding 

(60:0019838) (t, P<0.05) 

kinase inhibitor activity 

(GO:0019210) (t, P<0.05) 

peptidase inhibitor activity 

(GO:0030414) (t, P<0.05) 

protein kinase regulator activity 

(GO;00 19887) (t, P<0.05) 

transferase activity, transferring acyl groups 

(GO:0016746) (f, P<0.05) 

transferase activity, transferring glycosyl groups 
(GO:0016757) (J., P<0.05) 

monooxygenase activity 
(GO:0004497) (i, P<0.05) 

Cellular components 

extracellular matrix 
(GO:0031012) (t, P<10"*) 

intracellular 

(GO:0031012) (t, P<0.05) 

receptor complex (GO;0043235) (T, P<0.05) 



Biological processes 

response to organic substance 
(GO:0010033) (t, P<0.05) 

regulation of metabolic process 

(GO:00 19222) (t, P<10"^) 

regulation of catalytic activity 

(GO:0050790) (t, P<10"^) 

regulation of cellular process 

(GO;0050794) (t, P<10'^) 

regulation of cellular component organization (GO:0051128) 

(t, P<0.05) 

nucleobase, nucleoside, nucleotide and nucleic acid 
metabolic process 

(GO;0006139) (t, P<10"^) 

macromolecule biosynthetic process (GO:0009059) (t, P<10'^) 



gene expression 

(GO:0010467) (f, P<10"^) 

cellular biosynthetic process 

(GO:0044249) (t, P<10'^) 
programmed cell death 

(GO:0012501) (t, P<10"^) 

^ Specifically enricfied GO terms are fiighlighted in bold fonts 

The symbol for | represents the enrichment status of "enriched" or "depleted" of the GO term, respectively. 
' P values were Bonferroni-corrected for multiple tests 



transferase activity, transferring phosphorus-containing groups 
(GO:00 16772) (f, P<10"') 

hydrolase activity, acting on acid anhydrides (GO:0016817) 

(t, P<10''') 

protein domain specific binding 
(GO:0019904) (T, P<10^) 

substrate-specific transmembrane transporter activity 
(GO:0022891) (t, P<10"^) 

cytoskeletal protein binding 
(GO:0008092) (t, P<10"^) 

phospholipid binding 
(GO:0005543) (f, P<10"^) 

identical protein binding 
{GO:0042802) (f, P<10"^) 

receptor binding 

(GO:0005102) (f, P<10"'') 

heat shock protein binding 
(GO:0031072) (f, P<0.05) 

oxidoreductase activity, acting on peroxide as acceptor 
(GO:0015684) (t, P<10'^) 

extracellular matrix 

(GO:0031012) (f, P<10"^) 

extracelluar space 
(GO:0005615) (t, P<10"^) 

endomembrane system 
(GO:0012505) (t, P<10"^) 

membrane-bounded vesicle (GO:0031988) (t, P<0.05) 
cell surface 

{GO:0009986) (t, P<10'^) 
cell projection 
(GO:0042995) (t, P<0.05) 

midbody 

(GO:0030496) (f, P<0.05) 

response to organic substance 
(GO:0010033) (t, P<10"^) 

lipid metabolic process 
(GO:0006629) (f, P<10"'^) 

defense response 
{GO:0006952) (f, P<10"^) 

response to other organism 
(GO:0051707) (t, P<0.05) 

carbohydrate metabolic process 

(GO:0005975) (f, P<0.05) 

response to drug 
(GO:0042493) (t, P<10"^) 

cell-cell signaling 
(GO:0007267) (t, P<10"^) 

cytoskeleton organization 

(GO:0007010) (f, P<10'^) 

regulation of developmental process (GO:0050793) (t, P<10"'*) 



regulation of immune system process (GO:0002682) (t, P<10'^ 
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ExpDi^„ which represents the dissimilarity in expression 
profile, was calculated by 1 -Pearson's correlation coeffi- 
cient of the expression signals of the six tissues. ExpD^ao 
which represents the summed difference in expression 
levels, was calculated by the Euclidean distance of the 
expression signals of the six tissues. CTCF-binding sites 
identified by ChlP-seq experiments [25] in 13 non-can- 
cerous human cells (Table S5 in Additional file 1) were 
obtained from broadPeak data deposited as GSE30263 at 
NCBI GEO. We generated two sets of CTCF-binding 
sites. The number of overlapping CTCF-binding sites, 
which are CTCF-binding regions present in all examined 
cell types (Table S5 in Additional file 1), was determined 
using BEDTools [39]. The number of joint CTCF-binding 
sites, the union of CTCF-binding sites from all cell types, 
was determined using a custom Perl script. Enrichment 
analyses on GO terms were performed using FatiGO 
[40]. Partial correlation analyses were performed using 
modules of the "ppcor" package (v.1.0) [41] for R (http:// 
www.r-project.org/). 

Additional material 



Additional file 1: Supplementary tables S1-S5. 
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